Instability and Nonlinear Evolution of Narrow-Band Directional Ocean Waves 
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The instability and nonlinear evolution of directional ocean waves is investigated numerically 
by means of simulations of the governing kinetic equation for narrow-band surface waves. Our 
simulation results reveal the onset of the modulational instability for long-crested wave-trains, which 
agrees well with recent large-scale experiments in wave-basins, where it was found that narrower 
directional spectra leads to self-focusing of ocean waves and an enhanced probability of extreme 
events. We find that the modulational instability is nonlinearly saturated by a broadening of the 
wave-spectrum, which leads to the stabilization of the water-wave system. Applications of our 
results to other fields of physics, such as nonlinear optics and plasma physics are discussed. 
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Giant freak waves, or rogue waves, have been observed 
in mid-ocean and coastal waters [1|, in optical systems 
[5] , and in parametrically driven capillary waves |3] . The 
freak/rogue waves are short-lived phenomena appearing 
suddenly out of normal waves, and with a small probabil- 
ity [4j . The study of extreme gravity waves on the open 
ocean has important applications for the sea-faring and 
offshore oil industries, where they may lead to structural 
damage and injuries to personnel [T]. It is, therefore, 
very important to understand the physical mechanisms 
that lead to the formation of freak waves. Since the lin- 
ear theory cannot explain the number of extreme events 
that occur in the ocean and in optical systems, one has 
to account for nonlinear effects (e.g. wave-wave inter- 
actions) in combination with the wave dispersion. This 
can lead to the modulational instability (for water waves 
called the Benjamin- Feir instability O [6]), followed by 
focusing and amplification of the wave energy. 

Wind-driven waves on the ocean often have wide fre- 
quency spectra that are peaked in the direction of the 
wind [Tmi3j . The statistics of directional spectra for 
narrow-band gravity waves have also recently been stud- 
ied experimentally in water basins |14H16| . where it was 
found that sea states with narrow directional spectra 
(long-crested waves) were more likely to produce extreme 
waves. Examples of statistical models that govern col- 
lective interactions of groups of water waves are Hassel- 
mann's model |7j for random, homogeneously distributed 
waves and Alber's model [5j for narrow-banded wave 
trains. Wave-kinetic simulations in one spatial dimension 
have shown Landau damping and coherent structures 9J , 
and recurrence phenomena [10] for random water wave 
fields. In this Letter, we derive a nonlinear wave-kinetic 
(NLWK) equation for gravity waves in 2 + 2 dimensions 
(two spatial dimensions and two velocity dimensions) and 
carry out simulations to study the stability and nonlinear 
spatio-temporal evolution of narrow-band spectra waves 
that were observed in the recent experiments by Onorato 
and coworkers fT3]. The present NLWK model, which is 



similar to Alber's model is particularly suitable for 
studying the nonlinear dynamics of narrow-band water 
waves due to its relative simplicity. Similar nonlinear 
wave-kinetic equations also appear in the description of 
optical systems, photonic lattices, and plasmas jT7) . 

Deep water gravity waves are governed by the disper- 
sion relation ui = ^/gk, where g is the gravitational con- 
stant, k = y^fc^ + k'^ is the modulus of the wave vector 

k — kx^ + kyY, and x and y are the unit vectors in the 
X— and y— directions. Assuming surface displacements 
of the form 77 = (1/2) A(r, t) exp{—iujot + ikox)+ complex 
conjugate, where A is the slowly varying {\d/dt\ ^ utQ, 
|V| ^ ko) envelope, r — xx. + yy is the spatial coordi- 
nate, and ujq = y/gk^, the nonlinear interaction of water 
waves is governed by the nonlinear Schrodinger equation 
(NLSE) 

. (dA dA\ ^ d^A ^ d^A ^, ^,2 , „ 

where Vgr = duo/dkx = wo/2A;o is the group ve- 
locity, Dx = {ll2)d'^u:/dkl = -wo/8fcg and Dy = 
{l/2)d'^uj/dky = uJo/4:kQ are the group dispersion coef- 
ficients, and the nonlinear coupling coefficient is ^ = 
ujkQ/2. Introducing the two-dimensional Wigner trans- 
form [IH] 

/(r,v,t) = J ^*(R+,t)A(R_,Oe^^-(--''-^) d'X, 

(2) 

where we have denoted R± = r ± D • A and D • A = 
Dx\x^ + DyXyY, we obtain the evolution equation for 
the pseudo-distribution function / as 

X /(r, v', t)e'^-^^-^'U^v' cP\ = 0, 

where /(r, t) = J /(r, v, t) (Pv is the variance of the sur- 
face displacement (the wave intensity). The transfor- 
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mation ^ between ([I]) and ([s]) is valid in both direc- 
tions for a deterministic wave-train (corresponding to a 
"pure state" in quantum mechanics), with some restric- 
tions on the distribution function / [TB]; however, we 
are interested in the statistical properties of an ensem- 
ble of waves, and more general choices of / where the 
deterministic picture is abandoned [S]. In the absence of 
the nonlinear term in the left-hand side of (|3|, we have 
df/dt + V • V/ 0, which dictates that the wave en- 
ergy propagates in space with the group velocity v. Our 
model is valid for waves with v WgrX. The disper- 
sive properties of the wave are important for the non- 
linear wave-wave interactions between wave-packets that 
are modeled by the interaction integral in the last term 
in the left-hand side of ([S]). 

The velocity distribution can be related to the wave 
spectrum in frequency domain. Similar to Ref. |14| . we 
will use the model spectrum parameterized by the Joint 
North Sea Wave Project (JONSWAP) as [IT] 
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where ujp is the peak frequency, 7 is the peak enhance- 
ment parameter and ap is the Phillips parameter. Here 
7 is in the range 1-6 for ocean waves [14] , while ap is in 
the range 0.0081-0.1; the values 7 = 1 and ap = 0.0081 
gives the spectrum of fully developed wind seas [IH] , while 
the larger values are observed in water tank experiments. 
We will use ap « 0.025, 7 = 8 and a — 0.08, which are 
consistent with the Marintek water basin experiment in 
Refs. [MlllS]. Since the wave spectrum is concentrated 
around w = Wp, we will use ujq = ujp and kg = kp = uf^/g 
in the evaluation of and Dy in Eq. |3] 

The integral of the spectrum Q over all frequencies 
yields the variance of the surface elevation. While the 
variance of a monochromatic wave is |A|V2, from (§ we 
also have / f (Pv = \A\'^ /2. Hence, as initial conditions 
in our simulations, we will use / = /o(v) = Fq{v)G{9) 
where we have introduced polar coordinates = v cos(0) 
and Vy = V sin(0) in velocity space. We obtain Fq from 
the frequency spectrum Q by using the differential vari- 
ance dl = S{uj)duj — Fo(v) vdv, as 



F„{v) = S[oj{v) 



1 



dio 
dv 
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where we used that the group speed v of the wave 
packets is related to the wave frequency uj — \/gk via 
V = duj/dk = u}/2k = g/2a;, or uj(v) = g/2v. The 
directional spreading function is chosen as [T^ G{6) = 
Gocos^(6l/2) = Go[l -Hcos(6i)]^/V2^/2, where cos(6') = 



Vx/v, V = {vl+ w^)^/^, and Go 
stant \m such that /^^ G{e) dO - 



is a normalization con- 
: 1 . We note that G has 
a maximum at = and tends to a narrower distribution 
with an increase of the parameter N . 
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FIG. 1; (Color online) The time-evolution of the maximum 
intensity kll^ax for (a) N = 840 (black) (b) N = 200 (blue), 
(c) iV = 90 (red), (d) iV = 50 (green), (e) N = 2A (magenta), 
and the case of a narrow-band normal velocity distribution 
(the inset). The spatial distributions of wave intensity for 
(a)-(d) are shown in Fig. [s] at the times indicated here with 
arrows. 



35 
34 
33 
3.2 

3 1 






-0 




'0!p=0.03 








.^"0^=0.025 




e ' -0 ^> 


— 0- 


ap=0.02 



FIG. 2: (Color online) a) The linear growth rate ujj of the 
fastest growing wavemode and b) maximum kurtosis for = 
840 (black) N = 200 (blue), iV = 90 (red), iV = 50 (green), 
and A'^ = 24 (magenta), for ap = 0.02 (dashed line), ap = 
0.025 (solid hue) and ap ^ 0.03 (dash-dotted line). The solid 
line (ap = 0.025) corresponds to curves (a)-(e) in Fig. [l] 



Equation ([S]) can be cast into a numerically more con- 
venient form by employing the Fourier-transform in ve- 
locity space 



(6) 



/(r,r7,i) = 2 j f{r,^,t)e'^--d^v, 
which transforms Eq. ([s]) into 

^-*V^-V/+2i^[/(r+D.r,, t)-/(r-D.77, t)]/(r, q, t) = 0, 

where / = f{r,r],t)n=o/2. A similar equation was de- 
rived by Alber [5j, starting from the Davey-Stewartson 
equations for weakly nonlinear gravity waves. The nu- 
merical approximation of ([t]) is based on a method to 
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solve the Fourier transformed Vlasov equation [50]. Us- 
ing a pseudo-spectral method in space, the operator V is 
converted to multiplication by in, and the spatial shifts 
by ±D • f] in Eq. (ITj) are converted to multiplications by 
exp[±i(D • T]) ■ k], where k is the wave vector. The sys- 
tem was solved in a computational window moving with 
the group speed of the peak wave. We used a spatial 
domain of size x Ly = 100 x 500 resolved by 
Nx X A'j, = 32 X 32 intervals and with periodic boundary 
conditions, and a Fourier transformed velocity domain 
Lrjx X = 160ttv~,1 x IGOttv"^^ with Nrjx x N,jy = 80 x 80 
intervals, where Vph — ojp/kp is the phase speed of the 
peak wave. The velocity domain in our simulations is 

thus Vx^min — '^x — "^x^vaax and Vy^rnin — '^x — '^y.max 

where Vx,min = 0, Vx,max = '^'^ Nj^x / Lrjx = Sug^i and 
-Vy,min = Vy,max = T^Nny/Lr/y = Vgr. The simulation 
was initialized with the JONSWAP spectrum, where the 
Fourier integral ([6]) was evaluated numerically to obtain 
the spectrum in rj space. Random numbers of the order 
10~^ of the initial intensity was added to the solution 
in order to seed the modulational instability. The initial 
conditions give an intensity of / « 0.010 fc~^ uniformly 
distributed in space, which is compatible with the exper- 
iments of Onorato et al. [H]. To compare with the ex- 
perimental observations of Onorato et al. |14) , we carried 
out simulations for TV = 24, 50, 90, 200, and 840 corre- 
sponding to the Marintek experiment in Ref . p3] . They 
used ojp = 27r s~^ (1 Hz) and corresponding fcp = 4.1 m""'^, 
and a significant wave height Hg = 0.08 m, giving a wave 
intensity of / w 5 x 10^'' m^. 




0009 0.016 0.009 0.014 0.009 0.013 0.01 0.011 



FIG. 3: (Color online) The spatial distribution of the normal- 
ized wave intensity k^I for (a) A'' = 840 at f = 1.27 x 10^ Tp, 
(b) = 200 at t = 1.46 x 10^ Tp, (c) iV = 90 at t = 
1.81 X 10^ Tp, and (d) iV = 50 at f = 2.90 x 10^ Tp, corre- 
sponding to the curves (a)-(d) in Fig. [l] 

According to the analysis of Alber [8 , using a model 
two-dimensional normal spectrum, there are two condi- 
tions for the modulational instability: first, the modu- 
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FIG. 4: (Color online) The velocity distribution ojp/ of the 
wave energy, averaged over space, at i = (left column) and 
t = 3.2 X 10 Vp (right column), for (a) N = 840 (b) N = 200, 
(c) = 90, (d) = 50. Panel (e) shows the narrow-band 
normal velocity distribution aX t = Q (left) and t = 640 Tp 
(right). 



lational wavenumbers must lie within a certain direc- 
tional range (in Alber's case \Kx\ > \f2\Ky \ similar to the 
Benjamin- Feir instability), and second, the wave steep- 
ness (the wave amplitude multiplied by fcp) must be larger 
than the normalized (by the component of the spectral 
peak) spectral bandwidth. In our simulations, using di- 
rectional JONSWAP spectra, we observed the modula- 
tional instability and the self-focusing of the wave energy 
into localized wave packets for iV larger than 24. We 
measured the maximum value of the energy density in 
the simulation domain, and plotted its time evolution in 
Fig. [l] (the time is give in units of the peak wave period 
Tp = in/iOp). Initially, there is an exponential growth 
phase, reminiscent of the Benjamin-Feir instability for 
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monochromatic wave trains [5]. The modulational insta- 
bility is fastest growing for N = 840, and decreases with 
decreasing values of A^. For iV = 24 we do not observe 
any instability. For modulationally unstable cases, the 
exponential growth phase is followed by a nonlinear satu- 
ration of the instability, and finally a decrease of the max- 
imum energy density down to its initial background value 
/ ^ O.Olfc"^, as seen in curves (a)-(d) of Fig. [l] The inset 
shows a simulation with a narrow-band normal distribu- 
tion of the form / = 4w~^ exp[— 2(u^ -I- {v^ — Vgr)'^)/a^] 
with a = O.OA Vph, which yields the initial wave intensity 
/ = O.Olfc"^ that is similar as in curves (a)-(d). This 
case shows a rapidly growing instability to large ampli- 
tudes and then a decrease. The linear growth rate w/ of 
the instability for different values of N and a p was mea- 
sured from the data and plotted in Fig. [2ja). The growth 
rate is larger up to some limiting value for long-crested 
waves with N > 10^, while it approaches zero for smaller 



values of N. A growth rate of w/ = 1-2x10 



implies 



O.OOSwp, which is close 
,pUJp for monochromatic 



an amplitude doubling of the unstable wave in 50-100 
wave periods. The growth rate is sensitive to changes of 
ap and shows an increase/decrease of 50% with an in- 
crease/decrease of ap by 20%; this is consistent with a 
ratio of unity between the wave steepness and the spec- 
tral bandwidth, so that the system is weakly unstable. 
The strongly unstable case for the narrow normal distri- 
bution has a growth rate ujj 
to the limiting value [S] w/ = Ik^ 
waves. 

The kurtosis is traditionally [21] estimated by the for- 
mula A4 = 3 -|- 24fcpCr^, where a is the standard deviation 
of the surface elevation. (The factor 3 comes from the 
assumption of Gaussian statistics and the term 2Ak^a^ 
is a nonlinear correction to the Gaussian statistics.) As- 
suming that the wave field is ergodic, we have = (/), 
where (/) is the spatially averaged wave intensity. As 
noted in Ref. [IP , this formula underestimates the kurto- 
sis compared to the experimental values for narrow-band 
water waves, where an increase of the kurtosis was ob- 
served at later stages of the wave dynamics. Our model 
also conserves (/) and hence the formula predicts con- 
stant kurtosis. Taking into account that the wave-field 
is non-stationary and that the wave intensity varies in 
space (see Fig. 3), we, instead, estimate the kurtosis as 
A4 3(/^)/(/)^ + 24fcp(/), which assumes that the sur- 
face obeys Gaussian statistics locally everywhere. Using 
this estimate, we see in Fig. 2(b) that larger N gives 
larger kurtosis, in good agreement with experimental ob- 
servations [HHTB] . Figure [s] shows that the wave energy 
is concentrated into narrow bands, elongated along the y- 
dircction, which are propagating from left to right with 
speeds close to Vgr- At later stages, the wave-packets 
start to break up due to the two-dimensionality in space 



and the elongated bands of wave energy become more 
and more wiggled with the appearance of obliquely prop- 
agating waves, similar to those observed in Ref. [6]. For 
the modulationally unstable cases, the nonlinear interac- 
tion leads to a broadening of the distribution function in 
velocity space, as seen in Fig. [4j This, in turn, leads 
to a stabilization of the system via phase mixing of the 
wave envelopes [5, and a saturation and decrease of the 
maximum intensity shown in Fig. 1. 

To summarize, we have performed a series of kinetic 
simulations of narrow-banded water waves for different 
degrees of directional energy spectra. We observe an on- 
set of the modulational instability and self-focusing of 
the wave energy for waves with narrow directional spec- 
tra, leading to an increase of the estimated kurtosis. The 
modulational instability saturates via the occurrence of 
narrow wave-packets, which later disperse due to the 
broadening of the wave spectrum. Our simulation re- 
sults are in excellent agreement with observations from 
recent large-scale experiments in wave-basins |14H16j . 
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